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In bulk systems, molecules are routinely identified by their vibrational spectrum using Raman or infrared 
spectroscopy. In recent years, vibrational excitation lines have been observed in low-temperature conductance 
measurements on single molecule junctions and they can provide a similar means of identification. We present 
a method to efficiently calculate these excitation lines in weakly coupled, gateable single-molecule junctions, 
using a combination of ah initio density functional theory and rate equations. Our method takes transitions 
from excited to excited vibrational state into account by evaluating the Franck-Condon factors for an arbitrary 
number of vibrational quanta, and is therefore able to predict qualitatively different behaviour from calculations 
limited to transitions from ground state to excited vibrational state. We find that the vibrational spectrum is 
sensitive to the molecular contact geometry and the charge state, and that it is generally necessary to take more 
than one vibrational quantum into account. Quantitative comparison to previously reported measurements on 
pi-conjugated molecules reveals that our method is able to characterize the vibrational excitations and can be 
used to identify single molecules in a junction. The method is computationally feasible on commodity hardware. 



In recent years, vibrational excitations of single molecules 
have been measured in the inelastic tunneling regime with 
scanning tunneling microscopes (STMs),' in mechanical 
breakjunctions (MBJs)^'' and in the sequential tunneling 
regime (SET) with electromigrated breakjunctions (EBJs)."*"^ 
Measurements in EBJs on an oligophenylenevinylene deriva- 
tive by Osorio et alJ reveal a vibrational spectrum of 17 ex- 
citations in the sequential tunneling regime (see figure 5a). It 
has been shown that these modes are consistent with Raman 
(above 15 meV) and infrared (above 50 meV) spectroscopy 
data. However, the Raman and IR data show more peaks than 
are observed in the transport measurement. Moreover, the Ra- 
man and IR measurements were performed on polycrystalline 
samples and KBr pellets respectively, which do not reflect the 
conditions of the molecule in the junction. 

Theoretical investigations on vibrational excitations in the 
SET regime have so far mainly concentrated on small systems 
with only one vibrational mode*^ Chang et al}^ use den- 
sity functional theory (DFT) calculations on the C72 fullerene 
dimer to calculate all vibrational modes, restricting them- 
selves to transitions from the ground state in one charge state 
to the vibrational excited state of the other charge state. We 
have developed a computationally efficient method to calcu- 
late the vibrational spectrum of a sizeable molecule in the se- 
quential tunneling regime, based on first principles DFT cal- 
culations to obtain the vibrational modes in a three-terminal 
setup. This method takes the charge state and contact geome- 
try of the molecule into account and predicts the relative inten- 
sities of vibrational excitations. In addition, transitions from 
excited to excited vibrational state are accounted for by evalu- 
ating the Franck-Condon factors involving several vibrational 
quanta. Our method can therefore predict qualitatively dif- 
ferent behaviour compared to calculations that only include 
transitions from ground state to excited vibrational state. 

A schematic picture of a gateable electromigrated 
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breakjunction containing a single molecule is shown in fig- 
ure la. The couplings to the leads are given by Ti and Tr. In 
the weak coupling limit, Yi, TR,kT <sc A/s, Ec, the level spac- 
ing (Aii) and charging energy {Eq) of the molecule allow only 
one electron to tunnel onto the molecule at a time (sequential 
tunneling). The transport mechanism for this case is shown 
in figure lb. While tunneling on or off a molecule, an elec- 
tron can excite a vibrational mode, which may show up as a 
line running parallel to the diamond edges in the stability dia- 
gram (a plot of the conductance as a function of bias and gate 
voltage).'^ In this paper we calculate such stability diagrams 
with a rate equation approach and present a comparison with 
experimental results. 



Results and Discussion 

We have applied our rate equation method to three 
molecules with increasing length: benzenedithiol (see figures 
2 and 3), the oligophenylenevinylene derivative OPV-3 (see 
figure 4) and OPV-5 (see figure 5). All vibrational mode cal- 
culations are carried out using the Amsterdam Density Func- 
tional code""''^ with the Analytical Second Derivatives mod- 
ule.'** 

As a first example, we present our results for benzenedithiol 
adsorbed on gold. This example is used in order to demon- 
strate our method; In experiments, this system is generally 
not weakly coupled. We test the method by studying the in- 
fluence of the number of vibrational quanta, the charge state 
and the presence of gold contacts on the stability diagrams. 
The stability diagrams are calculated with a symmetric cou- 
pling to the leads of 1 meV, a bias (a) and gate (P) coupling 
of 0.5 and a temperature of 1 .6 K. The resulting stability dia- 
gram for the -1— >0 transition in bare benzenedithiol with one 
vibrational quantum (961 Franck-Condon factors) is shown in 
figure 2a. Of the 25 vibrational modes with energies below 
200 meV, eight excitation lines are visible belonging to the -1 
state, and seven belonging to the neutral state. For both states 
there are three excitation lines that do not continue all the way 
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FIG. 1 : (a) Schematic picture of a gateable junction containing a single molecule capacitively coupled to the left (Cl), right (Cr) and gate 
(Co) electrode, (b) The transport mechanism in the sequential tunneling regime. The distribution of the bias voltage over the leads is given by 



Cl+Cr+Cg 



and the gate coupling by fj ■■ 



Co 



Cl+Cr+Cg ■ 



to the diamond edge, but stop at the strong line at 43 meV. 

Taking two vibrational quanta into account (see figure 2b), 
reveals two new excitation lines for the neutral state, a strong 
line at 86 meV at a weak line at 180 meV (see the white ar- 
rows). Taking one more quantum into account (see figure 2c) 
adds a weak excitation line for the neutral charge state at 
129 meV, again indicated by a white arrow. For the 0— >+l 
transition (figure 2d), small changes are found. The arrows 
point to excitation lines that are absent in figure 2c. Com- 
pared to that diagram, the higher energy excitations (above 
100 meV) have shifted by several meV. 

In a junction, the molecule is bonded to the gold contacts. 
We have modeled this by adding two gold atoms on either 
side of the molecule. The resulting stability diagrams (with 
/ = 2) are shown in figure 3. These diagrams are quite differ- 
ent from those of the same charge state transitions in figure 2c 
and d. Depending on the charge state, five to eight excita- 
tions are visible below 75 meV, but no higher modes are ob- 
served. The electron-phonon couplings for the neutral charge 
state in the transition of figure 3a are shown in figure 7a. Two 
modes have a large electron-phonon coupling (with coupling 
strengths larger than 1), showing that it is necessary for this 
system to take more than one vibrational quantum into ac- 
count. 

The calculations show that only a few of the 30 vibrational 
modes of benzenedithiol are expected to be visible in trans- 
port measurements and that they are dependent on the charge 
state and sensitive to the contact geometry. For some modes in 
this molecule it is necessary to take more than one vibrational 



quantum into account. For example, the modes at 86 meV (for 
/ = 2) and 129 meV (for / = 3) are probably higher harmonics 
of the strong excitation at 43 meV. The fact that several other 
lines stop at this excitation shows that it is also necessary to 
take the Franck-Condon factors for excited vibrational state to 
excited vibrational state into account. 



The second molecule for which we have calculated the vi- 
brational spectiTim is OPV-3. As with benzenedithiol, the gold 
contacts are simulated by adding two gold atoms on either side 
of the molecule. The results for two charge state transitions 
with gold and one without are shown in figure 4. The calcula- 
tions take two vibrational quanta into account. Comparing the 
calculations to those on benzenedithiol indicates that OPV-3 
is less sensitive to the contact geometry. OPV-3 without gold 
has more modes at lower energies than benzenedithiol and the 
modes at higher energies are much less suppressed when the 
two gold atoms are added. Also, the electron-phonon cou- 
plings for OPV-3 are smaller than for benzenedithiol (see fig- 
ure 7b). These trends are not unexpected since OPV-3 is 
a larger molecule and the atoms will on average be further 
away from the leads, leading to a smaller sensitivity to the 
contact geometry. Also, since OPV-3 is conjugated, an extra 
electron will be delocalized over the entire molecule, and the 
atomic displacements will be smaller, resulting in a smaller 
electron-phonon coupling. In the case of OPV-3 we have per- 
formed several calculations with different contact geometries. 
We find that adding up to 19 gold atoms on either side of the 
molecule has no significant effect on the vibrational modes 



3 



tit'' 



a) 



b) 



I = 1 



1 = 2 



200 



200 





> 

E 




-200 




-200 



C) 



-200 




200 



d) 



1 = 3 



200 



S 




200 



^ 



1 = 3 




-200 



-200 



200 



-200 



-200 




200 



FIG. 2: Calculated stability diagrams of the -1— >0 (a-c) and 0— »-l-l (d) transitions in benzenedithiol with increasing number of vibrational 
quanta. The arrows point to the main differences between the diagrams (see text). Since the calculation is symmetric in the bias voltage, they 
are only shown for positive bias. 



above 20 meV. 

The calculated stability diagram of the third molecule, 
OPV-5, is shown in figure 5b. The temperature and coupling 



' We have also performed measurements on vibrational excitations in OPV- 
3. However, broadening of the lines due to large couplings to the leads pre- 
vents us from obtaining measurements with sufficient resolution to make a 
quantitative comparison to the calculations possible. The measurements do 
show the same trends as the calculations. None of the samples show any 
excitations above 30 meV, and only a few below. 



parameters of the calculation are chosen to be the same as in 
the experiment (figure 5a). Although we are unable to deter- 
mine the charge states in the measurement, the fact that the 
degeneracy point is the first at a negative gate voltage sug- 
gests a -1— »0 transition (see also figure 3 in^). In the calcula- 
tion, the non-conjugated dodecane sidearms of the measured 
molecule are omitted. These arms are not expected to influ- 
ence the electronic transport and will most likely only affect 
the low-energy vibrational modes. As with benzenedithiol, 
the contacts are modelled by adding two gold atoms on either 
side of the molecule. The calculation takes one vibrational 
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FIG. 3: Calculated stability diagrams of the -1— »0 (a) and 0— > + 1 (b) transition in benzenedithiol with two gold atoms on either side to simulate 
the leads. Two vibrational quanta are taken into account. 
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FIG. 4: Calculated stability diagrams for OPV-3 with (a and b) and without (c) gold for two charge state transitions. The calculations take two 
vibrational quanta into account. 



quantum into account. 

Figure 6 shows a t/Z/t/V-trace along the diamond edge of 
the neutral charge state in figure 5b. However, since this is a 
smaller calculation, three vibrational quanta can be taken into 
account. The peaks in this figure correspond to the excitation 
lines in the calculated stability diagram. In the experimental 
stability diagram, a background conductance makes it diffi- 
cult to resolve all excitation lines at the same color scale, but 
close inspection reveals 17 modes in the energy range below 
125 meV (see figure 6 and table I in^). The energies of the 
excitations in the measurement are determined from the bias 
voltage at which they cross the diamond edge. Broadening 



due to the temperature and the leads introduces an uncertainty, 
indicated by the horizontal bars in figure 6. 

Figure 6 reveals a close match between the experiment and 
the calculation for the modes between 10 and 80 meV. The 
calculation shows several small peaks in this range not ob- 
served in the measurement. It should be noted that in the rate 
equation approach, broadening of excitation lines is solely due 
to temperature. Broadening due to the couplings to the leads 
is not accounted for. Calculations which do take this broaden- 
ing into account show that these small peaks are smeared out 
and the calculation and measurement show the same number 
of peaks in the aforementioned range. 
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FIG. 5: (a) Measured stability diagram of a junction containing OPV-5.' (b) Calculated stability diagram of an OPV-5 molecule with one 
vibrational quantum, (c) The configuration of the OPV-5 molecule in the calculations. The dodecane sidearms of the measured molecule (see 
figure la in') are omitted and two gold atoms are added on either side to simulate the leads. 



The inset in figure 6 shows the same calculation, but with 
the gold atoms omitted. Comparison with the measurement 
shows a large discrepancy for excitations below 50 meV. It is 
clear from this figure that the addition of two gold atoms on ei- 
ther side of the molecule can already account for most of the 
influence of the contact geometry on modes above 10 meV. 
The charging energy of an OPV-5 molecule in a junction is 
an order of magnitude smaller than the difference between the 
ionization energy and electron affinity of the molecule in the 
gas phase, ''^ probably due to screening in the leads. While 
this effect changes the energies of the orbitals, the shape of 
the orbitals, and therefore the electron density will remain rel- 
atively unaffected. Since the Franck-Condon factors primarily 
depend on the difference in electron density between different 
charge states, we have chosen not to take image charges into 
account in the calculations. 

While this effect changes the energies of the orbitals, the 
shape of the orbitals, and therefore the electron density, on 
which the Franck-Condon factors primarily depend, will re- 
main relatively unaffected, and we have chosen not to take 
image charges into account in the calculations. 

The omission of the sidearms in the calculation lowers the 
mass of the molecule, which might explain the discrepancy 
between the calculation and the measurement for the modes 
below 10 meV, which involve motions of the whole molecule. 
Also, the contact geometry in the measurement is unknown, so 
any mode involving a significant distortion of the gold-sulfur 



bond is expected to be inaccurate. Like OPV-3, the vibra- 
tional spectrum of OPV-5 is less sensitive to the charge state 
and contact geometry than benzenedithiol and the electron- 
phonon couplings are smaller (see figure 7c). As in the case 
of benzenedithiol and OPV-3, the calculation of OPV-5 pre- 
dicts the intensity of the excitation lines to be much weaker 
above 30 meV than below. This is also observed in the mea- 
surement. The intensities gradually increase for energies up to 
30 meV, after which they suddenly drop, a trend also visible in 
the electron-phonon couplings. For excitations above 80 meV, 
the low intensities make a quantative comparison between the 
measurement and the calculation difficult. 

Most of the vibrational modes have electron-phonon cou- 
plings below 0. 1 and are not expected to give rise to extra exci- 
tation lines when another vibrational quantum is taken into ac- 
count. The modes at 17 and 27 meV, with coupling strengths 
of respectively 0.6 and 0.7, are expected to give rise to exci- 
tation lines at 34, 51-54 and possiby 81 meV. These lines are 
indeed observed in the measurement and the calculation (see 
figure 6). 

It should be emphasized that in figure 6, all visible vibra- 
tional excitations, for both the calculation and the measure- 
ment, are shown. Comparing the spectrum to Raman and IR 
spectroscopy data reveals a close match,^ but the optical spec- 
tra predict many more modes not observed in the measure- 
ment and calculation. The calculation predicts only a handful 
visible excitations out of a total of a 129 vibrational modes 
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FIG. 6: Calculated dl /dV-trsLce of the diamond edge of the neutral charge state in figure 5b, taking three vibrational quanta into account. The 
inset shows the same calculations, but with the gold atoms omitted. All measured excitations in this energy range (see figure 5a) are shown. 
The uncertainties in the measured energies are indicated by the horizontal bars. 
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FIG. 7: Dimensionless electron-phonon couplings (A) for the vibrational modes of the neutral charge state at the -1 — » transition for (a) 
benzenedithiol, (b) OPV-3 and (c) OPV-5. All calculations include two gold atoms on either side of the molecule to simulate the leads. Analysis 
of the atomic displacements shows that primarily modes that distort the n - n overlap give rise to a non-zero electron-phonon coupling. 



under 150 meV. Our method is thus able to provide what 
we might call 'selection rules' for vibrational excitations in 
single-molecule junctions. 



Conclusions 



In conclusion, our results show that the vibrational spec- 
trum of a single-molecule in a junction is sensitive to the 
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contact geometry and charge state, although this influence be- 
comes smaUer for larger molecules. Contrary to Raman and 
IR spectroscopy, calculations can take these influences into 
account, provide 'selection rules' and predict the relative in- 
tensity of excitation lines in transport measurements. Our cal- 
culations also show that it is necessary to take more than one 
vibrational quantum into account for small molecules, but that 
due to decreasing electron-phonon couplings this becomes 
less important for larger molecules. Finally, our method is 
computationally efficient. All Franck-Condon and transport 
calculations have been performed on an HP xw9300 work- 
station, with the largest calculation (over 54 million Franck- 
Condon factors in OPV-3 for 250.000 bias and gate points) 
taking just over four hours. 



Methods 

Rate equations. Our method to calculate stability diagrams 
is based on the rate equation approach, which is generally 
used in the sequential tunneling regime. The central quan- 
tity in this approach is the vector of occupation probabilities 
Pn,<T,v for states with charge quantum number n, spin quantum 



number cr and vibrational quantum number v. These proba- 
bilities change by transitions from one state to the other with 
rates R„^a-,v^n',a-'y- The time evolution of the occupation prob- 
abilities as a function of these rates is described by the master 
equation: 



dt 



(1) 



n' ,cr' y ,tn,(7,v 



The rates are given by Fermi's golden rule: 



R 



1 1/ , , ,1 \|2 
',o-',v' ~ ^ |\" ^ |//| nerval p„'_ 



(2) 



These rates contain contributions from the electronic and nu- 
clear wavefunctions. The electronic contributions are de- 
scribed by the couplings to the leads and the Fermi function, 
and the nuclear contributions are described by the Franck- 
Condon factors, which we will discuss below. 

At low bias, only two charge states are relevant: the initial 
state (n,(T',v'), where the level in the bias window (see fig- 
ure lb) is unoccupied, and the final state (« +l,cr, v), where it 
is occupied. The rates for these states are:'-''^^ 



R 



«,o"',v'— >/I+l,0", 



(3) 



Kyy^n+i,a,v =Fy,y-^f[E„+i,r^y - E„^^,y - epVg -e(a-l) Vt) 

r 



17+1, tr,v—>n,a"',v' 



^Py'y-T [l -f{En+i,a,v - En,o-'y - epVg - e (1 - a) V/,)] 



r 



where Fy'y are the Franck-Condon factors and / is the Fermi 
function. E„+\^a-,v-E„^,r'y is the energy difiference between the 
initial and final state. This difference is composed of the level 
spacing, the charging energy, the vibrational reorganization 
energy and the vibrational energy of the states v' and v. For 
a single degeneracy point, all but the latter of these terms can 
neglected by choosing a suitable reference point for Vg. 

With the rate equations in place, we write the master equa- 
tion ( 1 ) in matrix- vector form: 



dt 



MP 



(4) 



where P has elements P„^,r'y ™d P„+i^a-.v M is the 2N x 
2N rate matrix, where numbers the vibrational states. Its 



elements are given by 

-.iV 

-1 

R 



Mij 



Zjv=1 ^n.o"',/— >n+l,iT,v 



Zjv' = 1 '^n+\,o-J-N^n,cr' y 
^n+\ ,crJ-N^n,o-' J 
t^n.cr' j'— >/!+ 1 ,o-,i-N 





if / - j and /, j < N, 
if / - j and /, j > N, 
ifi<N and j > N, 
ifi>N and j < N, 
otherwise. 



(5) 



To calculate the current we need the stationary occupation 
probabilities, i.e. ^ = 0, which can be obtained by calculat- 
ing the null space of M, with the condition that all elements of 
P are non-negative and Tjn.o-.v Pn.o-.v - 1- Once the stationary 
occupation probabilities and the rates are known, the current 
can be calculated by summing over the total rate through one 
of the leads. For the left rate this becomes:**'^' 

N N 

t - ^ ^ ^ {P n.cr' yRfi^a-' y ^n+\,a-,v ~ ^"+l.o".''^«+l,o-,v^n,o-',v') 

v' = \ v=\ 

(6) 

Franck-Condon factors. The Franck-Condon factors 
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(Fyiy) in equation (3) are a measure of the probability that 
a tunneHng event will be accompanied by a vibrational ex- 
citation. When an electron tunnels on or off a molecule, 
the change in electron density will shift the equilibrium po- 
sition of the nuclei and possibly cause a transition to a dif- 
ferent vibrational excited state. The probability for this to 
happen is equal to the square of the overlap integral of the 
vibrational wavefunctions in both charge states. To calcu- 
late the overlap integral, the normal coordinates of one charge 
state have to be expressed in those of the others. This pro- 
cedure is known as the Duschinsky transformation.^"*'^^ This 
transformation yields the Duschinsky rotation matrix and a 
mass-weighted displacement vector fc. The latter can then 
be used to calculate the dimensionless electron-phonon cou- 
plings Ai = ^A:,-.22.26 

Two methods are mainly used for calculating the Franck- 
Condon factors from the Duschinsky transformation and the 
frequencies. One is the generating function method of Sharp 
and Rosenstock^^; the other is the recursion-relation method 
of Doktorov et al?^ We have used the latter as implemented in 
the two-dimensional array method of Ruhoff and Ratner.^^-^" 
This method returns an NxN array of Franck-Condon factors, 
where = is the number of permutations of up to / 
vibrational quanta over a vibrational modes. 

Relaxation rates. Since Franck-Condon factors represent 
probabilities, the elements of each row and each column of 
this matrix add up to 1 for I — > oo. Typically, for a large sys- 
tem with many Franck-Condon factors there are many factors 
where Fy>y <k 1 . This can present numerical difficulties when 
calculating the stationary occupation probabilities. The rates 
are proportional to the Franck-Condon factors and if all the 
rates into and out of a certain level are (very nearly) zero, any 
occupation of that level will be stationary, resulting in an in- 
finite number of solutions for the stationary occupation prob- 
abilities. There are two ways to prevent this from happening. 
The first is to take intramolecular vibrational relaxation into 
account. We have implemented a simple relaxation model in 
which a single relaxation time r is used for all states: 



dt 



where 



Z 



P f > >R f f f - P R > f > 

-* n\(T' y ^^n' ,(r'v -^n.cr.v ^ n.cr.v^^n.cr.v^n' .a' y 



n' ,cr' y i^n,o-,v 

1 ' 



p _ p"'' V p 



r^eq 



(7) 



(8) 



is the equilibrium population according to the Boltzmann dis- 
tribution. This term can be included in the rate matrix by 
adding the relaxation matrix: 



'J 



1 



peq 

peg 

n,cr' J 
pel 

'^n+\,cr,i-N 
peg 





if / = /' and /, j < N, 
if / ^ j and /, j < N, 
if ; - j and /, j > N, 
if / j and /, j > N, 
otherwise. 



(9) 



For sufficiently small relaxation times, the previously men- 
tioned stationary states will decay to the ground state and there 
will be only one solution. 

Iterative solution. The second approach is to calculate the 
null space iteratively by starting from the equilibrium popula- 
tion. Since both the equilibrium population of higher energy 
states and the rates to those states are zero, they will never 
become populated and the method will converge to the phys- 
ically correct solution. This approach has several additional 
benefits. The rate matrix for most realistic systems will be 
sparse and an iterative method can make use of this and scale 
better than a direct method. Also, the stationary population at 
neighbouring bias and gate points will be the same unless a 
new state becomes available, so by using the previous popula- 
tion as a starting point, most points will only need a single it- 
eration to converge. We have implemented a direct method us- 
ing singular value decomposition and an iterative method us- 
ing a Jacobi preconditioned biconjugate gradient method. The 
second implementation is generally several orders of magni- 
tude faster It turns out that implementing a combination of 
both approaches yields a one-dimensional null space of the 
rate matrix, even for larger molecules with several vibrational 
quanta. 
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